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PACS 87. 10. Hk - Lattice models in biological physics 

PACS 87. 16. Ka - Filaments in subcellular structure and processes 

PACS 87 . 16 . Qp - Flagella 

Abstract - The growth of bacterial flagellar filaments is a self-assembly process where flagellin 
molecules are transported through the narrow core of the flagellum and are added at the distal 
end. To model this situation, we generalize a growth process based on the TASEP model by 
allowing particles to move both forward and backward on the lattice. The bias in the forward 
and backward jump rates determines the lattice tip speed, which we analyze and also compare 
to simulations. For positive bias, the system is in a non-equilibrium steady state and exhibits 
boundary-induced phase transitions. The tip speed is constant. In the no-bias case we find that 
the length of the lattice grows as N(t) oc yi; whereas for negative drift N(t) oc hit. The latter 
result agrees with experimental data of bacterial flagellar growth. 
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£>~» Bacterial flagella act as motility devices that allow bac- 
i-Jh teria such as Escherichia coli and Salmonella to swim and 
i ^*t o respond to chemical stimuli by performing chemotaxis 
|T]. A flagellum mainly consists of a long helical hollow 
structure with a length of up to 20/ito and a typical outer 
J^l diameter of 0.02/Km. The growth mechanism of the flagel- 
lar filament is a self-assembly process. Flagellin molecules 
are transported through the hollow core of the filament 
CNj and are added one by one at the tip of the filament. An 
important quantity of such a growth process is the time 
i-H dependence of the filament length N(t). A recent model 
CN|that treats flagellin molecules as diffusing particles in a 
r j"^ single-file process established a growth function N(t) oc \ft 
[2J. However, the only available experimental data to our 
knowledge show a logarithmic growth N(t) oc hit [3J. To 
model the growth process and to try to verify the export- 
ed mental data, we make use of a variant of the well known to- 
tally asymmetric simple exclusion process (TASEP). This 
model takes into account the important fact that flagellar 
growth happens in a single-file process meaning flagellin 
molecules cannot pass each other. 

The TASEP itself has developed into a paradigmatic 
model of non-equilibrium physics - In particular, 
TASEP models with open boundaries showed to be very 
fruitful for the study of non-equilibrium steady states, 
i.e., states that are characterized by non- vanishing cur- 
rents. These systems have been used to model the mo- 
tion of ribosome along mRNA [TO] , molecular motors along 
microtubule filaments or that of cars on a highway 
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Fig. 1: Growing asymmetric simple exclusion process with open 
boundary conditions. The labels indicate rates at which par- 
ticle hops can occur. At site 1 particles transform into a new 
lattice site with rate 7, which then becomes the new site 1. 



(T2ll3| . just to name a few applications. All of these mod- 
els consider lattices with a fixed number of lattice sites 
N. However, there are recent approaches that general- 
ize the TASEP to dynamically extending exclusion pro- 
cesses. One approach considers the TASEP with a fluc- 
tuating boundary that can be pushed away by particles 
on the lattice Q3] , another one takes into account a death 
probability of the leading particle [TOJ[TO]. In a different 
variant of the TASEP, particles that reach the end of the 
lattice can transform into new lattice sites resulting in a 
dynamically extending lattice [T7]. Initially, the extend- 
ing TASEP was formulated to model fungal hyphal growth 
|18j . In the reminder of this paper we will generalize this 
model to a partially asymmetric model where particles can 
move both forward and backward on the lattice. In par- 
ticular, we calculate the growth function N(t) and discuss 
implications for the bacterial flagellar growth. 

The model is defined in the following way (see fig. I for 
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a schematic) . Particles enter the lattice with rate a at the 
boundary on the right. They either move to the left with 
rate q or to the right with rate p so that q+p = 1. Finally, 
at the tip of the lattice at site 1 they transform into a 
new lattice site with rate 7. A crucial feature of TASEP 
models is that particles can only move if the target site is 
empty. This hard-core interaction results in a nonlinear 
relationship between current J and density p, which itself 
is the origin of the many interesting features of the TASEP. 

In a first step we study systems where particle hops 
possess a positive drift q > p. To begin the analysis, we 
write down the averaged rate equations for occupancy n% G 
{0, 1} of site i: 

(ni) = q((l -ni)n 2 ) -n 2 )n 1 ) -7(711) , 

(n 2 ) = q((l - n 2 )n 3 ) - q((l - m)n 2 ) - - n 3 )n 2 ) 

+p({l - n 2 )ni) - j(nin 2 ) , 
(rij) = q{{l - ni)n i+ i) - q{{! - Ui-ijUi) 

-p((l - n. l+1 )n.i) +p((l - m)ni-i) i > 3 . 

+7(ni(rij_i - rij)) , 

The term 7(^1) in the first equation describes the growth 
of the lattice. Further terms proportional to 7 occur since 
after each growth event the lattice sites are renumbered. 
This corresponds to working in the frame of reference of 
the tip so that the new site at the left boundary becomes 
site 1. To tackle this system of coupled equations, we 
employ a mean-field approximation by setting (niUj) = 
(rii}(nj) = piPj, where pi is the density of site i. We 
now identify the particle current Jjj from site i to its 
neighboring site j as 

Ji,o = IPi , 

>h.\ = q{i - pi)p2 -p(l - P2)P\ , (l) 

Ji+i,% = 9(1 - Pi)Pi+i - P(l - Pi+i)(>i - IP1P1 , * > 2 . 

In a steady state, the current J is constant throughout 
the lattice. Moreover, since particles at the boundary on 
the left are used solely for the growth of the lattice, the 
current must exactly match the tip speed J = v = jpi. 
This allows us to solve the set of equations ((TJ) for the 
density throughout the lattice 

Pi-', P2- q ^^ [q _ pV 

(2) 



The last equation shows two fixed points (see fig. [2]) upon 
setting pi = Pi+v 

q-p-v fq-p-v\ 2 v 

where p_ is stable and p+ is unstable as indicated in fig. 
HJa). Iterating eq. (0 results in four possible density pro- 
files [sec fig. [UJb)]. Whereas the low-density profiles relax 




Pi 1 i N 



Fig. 2: Sketch of possible solutions to recurrence relation ([2} 
on the left and corresponding density profiles on the right. p_ 
is a stable fixed point whereas p+ is unstable. 

towards the stable fixed point p_, the high-density pro- 
files start in close vicinity of the unstable fixed point p+ 
and deviate from it when the other end of the lattice is 
reached. Upon changing bias q and boundary conditions 
a, 7, it is now possible to identify three different phases. 
For low input rate a at lattice site N, the system is in the 
low density (LD) phase, whereas low growth or release rate 
7 at lattice site 1 leads to the high density (HD) phase. 
The crossover between low density and high density phase 
is governed by a first-order transition, where both phases 
coexist. The locations of the phases as determined by 
the mean-field theory are similar to the refined a, 7 phase 
diagram in fig. [31 which we will discuss below. When 
both rates a and 7 are high, the fixed points of eq. ^ do 
no longer exist and the system is in the maximal current 
(MC) phase. The density profile is then no longer con- 
trolled by the boundary conditions a and 7. Instead, the 
bottleneck in the system is given by the bulk dynamics, 
namely due to the implicit hard-core interaction, which 
restricts the current or tip speed to its maximum value 
Jmc — vmc — (3 — 2^/2){q — p). It is determined from 
cq. ([3]) by setting the root to zero, i.e., when the two fixed 
points merge. Decreasing the hopping rate bias q below 
1 slows down the bulk dynamics which results in an ex- 
panding maximal current phase in the phase diagram. 

In order to determine the phases just discussed within 
the formulated mean-field theory, we extended and gener- 
alized the procedure of ref. [T7] to q < 1. We shortly sum- 
marize the reasoning. Setting pn-i = Pn m the mean- 
field rate equation for the boundary on the right, 

Pn = p(l - Pn)pn-i - <7(1 - Pn-i)pn + a(l - p N ) , (4) 

yields the boundary condition piq = a/ \q — p). In the 
LD phase, pm equals the stable fixed point p- of eq. ^ 
which gives an expression for the LD current or tip speed 
Jld = Vld- However, the iteration of eq. @ starting at 
p 2 only converges to the stable fixed point p_ if it starts 
below the unstable fixed point p + , i.e., if p 2 < p + . Using 
the expression for J^jj = vld, this condition leads to a 
phase boundary against the HD phase. In the HD phase 
the tip velocity vhd and the high density value p+ follow 
by setting p 2 = p + . From p 2 = p + the density then has 
to decay towards pn = &/(q — p) &t the boundary on the 
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right. This sets up a stationary shock profile at the phase 
transition between LD and HD phase. Close to this phase 
transition in the HD phase, the shock moves away from 
the tip but never reaches the other boundary of the lattice. 
Only further away from the phase transition, docs the high 
density expand throughout the lattice except in a narrow 
region close to the boundary on the right. 

Monte Carlo simulations of our model were done using 
a random sequential update where sites are picked ran- 
domly one after another and updated according to rates 
a, 7, q. One time step consists of updating all N lattice 
sites. The tip velocity v is then determined as a func- 
tion of either a or 7 for a given value of q. We localize the 
HD-MC or LD-MC phase transition when v reaches a con- 
stant maximum value which indicates the MC phase. As 
already indicated, the resulting phase diagram obtained 
by the mean-field approach does not coincide very well 
with results from Monte-Carlo simulations. An improved 
mean-field calculation that takes into account correlations 
between site 1 and 2 by setting {n\n2iii) = (n\n%)ni, pre- 
dicts a more precise phase diagram. This was already ob- 
served for 5=1, i.e., in the totally asymmetric case [17] . 
In contrast to the mean-field approach discussed above, 
iteration of eq. ([2} then starts at site 3 instead of site 
2. In particular, to eliminate (711712) in the equation for 
P3, one needs to evaluate the rate equation for the state 
(712(1 — ^l)}- The resulting phase diagram is shown in 
fig. [3] and compared to Monte-Carlo simulations. The im- 
proved mean-field calculations correctly predict the phase 
boundaries between LD and MC phases for various q val- 
ues whereas the boundaries between HD and MC phases 
still show some deviations from simulations. This was al- 
ready observed in [17]. Since determining the HD-MC 
boundary involves the starting value of iteration @, the 
result depends strongly on the number of site correlations 
that are taken into account at the tip. Figures Eta) and 
(b) show density profiles at constant a, 7 that convert 
from the respective LD or HD phase at q = 1 to the MC 
phase at q — 0.8 and 0.6. 

For our purposes the main conclusion is that the system 
reaches a non-equilibrium steady state characterized by a 
constant current or tip speed v. Hence, the lattice grows 
linearly in time: N(t) = vt. This is valid for all three 
phases and the actual value of v depends on the parame- 
ters a, 7 and q. With v M c = (3 - 2\/2)(q - p) in the MC 
phase, we obtain 



N(t) = (3-2V2)(q-p)t, q>p 



(5) 



which agrees with simulations [sec fig. Eta)]- Note that 
Eq. (O predicts N(t) —> for decreasing bias q —> 0.5. 
However, Monte-Carlo simulations show a non-vanishing 
growth of the lattice for q = p. So, we have to study this 
case separately. 

When q = p = 0.5, the bulk of the lattice is not driven 
anymore, i.e., the condition of detailed balance is fulfilled 
in the bulk but not at the boundaries. The system is still 




Fig. 3: Phase diagram of extending TASEP. Monte Carlo sim- 
ulation (symbols) and improved mean-field calculations (lines) 
are shown for q = 1.0, 0.8, and 0.6. Horizontal lines are boun- 
daries between HD and MC phases. Vertical lines separate LD 
from MC phases. Bullets: parameters of density plots in figs. 
Eta), (b). 



kept out of equilibrium. For q = p, eqs. © reduce to 
p 1 =v/'y, p 2 = 2v + v/~f, 

p l+1 = Pi (l + 2v) + 2v , i>2. 



(6) 



The recurrence relation does not have a fixed point with 
p > and the density grows to infinity, lim pi = 00. To 

i— >oo 

bypass this behavior, we impose that the density reaches 
the physical maximum at the boundary on the right, pjy = 
1. For N ^> 1 and with p^ = 1 one readily shows that 
pi = w/7 -C 1 and therefore negligible in eqs. Then, 
the recurrence relation of (|6|) is solved by 



Pi = (1 + 2V) 1 - 1 - 1 



(7) 

The condition p^ = 1 gives 2 = (1 + 2v) Nt - t \ which in- 
dicates that v is time-dependent now. By identifying v(t) 
as dN(t)/dt, we can solve for N(t) and arrive at 



N(t) = \A^2~Vi 



for 



t > 1 



(8) 



This result is confirmed by simulations as indicated in fig. 
Eta). It is independent of a and 7, i.e., the MC phase fills 
the whole phase diagram. It also shows that the current 
in the lattice is ohmic, J = v = dN(t)/dt oc 1/N(t). 

We note that the growth function of eq. (0) is consis- 
tent with a continuum approach for this particular system. 
The non-extending TASEP with fixed N, open bound- 
aries, and q = p reduces to the diffusion equation in the 
continuum limit. Hence, to obtain a continuous model of 
the extending TASEP with q = p, one has to formulate 
the diffusion equation on a growing domain. Similar to 
writing a reaction-diffusion equation on a growing domain 
□SI, we find M 



dp 
dt 



1 d 2 p 1 dN(t) 



2N(t) 2 dx 2 N(t) dt 



€[0,1], (9) 
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Fig. 4: Monte Carlo simulations of density profiles for different 
phases: (a) LD phase with crossover to MC phase (q — 0.8, 0.6), 
(b) HD phase with crossover to MC phase (q = 0.8, 0.6) [Bullets 
in fig. |3]give the location in the phase diagram], (c) q — p, and 
(d) q < p. 



where we replaced pi(t) by p{x, t) and the growing domain 
is mapped on the interval [0, 1]. This equation only has a 
stationary solution in the case of diffusive growth N(t) oc 
\ft. For dp/dt = and using growth function (O, we 
arrive at the time-independent equation 



J_fp 
ln2dx 2 



-p = 0. 



(10) 



Under the assumption of Dirichlct boundary conditions 
p(0) — and p(l) = 1, where the tip of the lattice is 
located at x = 0, the stationary solution reads 



p{x) = 



sinh • 



(11) 



Monte Carlo simulations confirm this density profile (see 
fig. |U(c)). The deviation from a linear density profile, i.e., 
the solution to the non-extending TASEP with q = p, is 
rather small. 

When q < p, the particles in the bulk prefer to hop to 
the right, even though the boundary conditions enforce a 
current from right to left. As a result, the growth func- 
tion ([5]) becomes negative indicating that the mean-field 
calculation is not applicable in this case. Studies of the 
non-extending TASEP with q < p exist [21][22] and exact 
solutions revealed a reverse bias phase where the parti- 
cle current decreases exponentially with system size N. 
In our case, we expect that the lattice grows very slowly 
for q < p. Furthermore, most of the particles will stay 
in a domain at the right boundary, whereas near the left 
boundary most of the sites are vacant. Monte Carlo sim- 
ulations confirm that there is a domain interface between 
an empty domain and a "jammed" domain. As fig. HJd) 
shows, this domain wall is roughly located in the middle 
of the lattice. Now, whenever the left-most particle in the 




Fig. 5: Trajectory of a biased random walker with fixed wall 
at the origin. The random walker starts at x — 0. The height 
of an excursion is called di, the length is called n. 

lattice reaches the tip, the lattice will extend by one site. 
After such a growth event, the next particle in the lattice 
becomes the left-most particle. Therefore, the growth of 
the lattice is equivalent to finding the probability that a 
biased random walker reaches a maximum distance from 
a bounding wall. In this picture, the wall is given by the 
domain interface between the empty domain on the left 
and the "jammed" domain on the right and the maximum 
distance is the distance of the lattice tip from the domain 
interface. Figure [5] depicts a trajectory of a biased random 
walker with fixed wall at the origin. The random walker 
is on excursions that last for a timespan of tj and have a 
maximum distance di from the origin. Ref. |23j determines 
the probability for di to be smaller than some threshold c 
in the long-time limit as 



P(d z < c) = 1 - e- c / K , 



(12) 



where k characterizes the step-length distribution W{X{) 
of the random walker by the condition 

(e x ' /K ) = 1 . (13) 

In our case, the respective probabilities for hops to the 
left or right, i.e., away from or towards the origin, are 
W(Xi = +1) = q and W(X, = -1) = p. The condition 
(fT3)l gives 



qe 1/K +pe- 1/K = 1 , 



-1/k 



\vh irti results in 



In? 



(14) 



At time t the random walker made t/fc) excursions and 
the maximum M t is given by 



M t =max(di,d 2) ...,d t /( n )) 
The probability that M t < c becomes 
P(M t < c) 



(15) 



P{d l < c)-P{d 2 < c)...P(d t/{Ti) < c) 

t/(n) 



= 1-e 



Replacing c by c + k In t yields 



P(M t < c + Klni) = 1 



-c/k 



t/{n) 



(16) 
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Fig. 6: Lattice growth determined from Monte Carlo simula- 
tions, (a) Upper 3 curves: positive drift with q > p, bottom 
curve: no drift with q = p — 0.5. (b) Negative drift with q < p. 



and in the limit of large t/(ri), one obtains the Gumbel 
distribution [Ml 



P{M t <c + K\nt)=exp(-e- c / K /(n)) 



(17) 



This means that the probability for a maximum excursion 
stays constant in time, when in the long-time limit the 
maximum grows as 



M t = K.lnt + c 



(18) 



where c is an undetermined constant. Since the domain 
interface is situated approximately in the middle of the 
lattice, the growth function of the whole lattice is then 



In - ) 



Q < P , 



with tip speed 



dN(t) N 



(19) 



(20) 



Interestingly, the dependence of v on N agrees quali- 
tatively with the dependence of the current on N in 
the non-extending TASEP with negative drift [2lH2"2] . 
Simulations of an extending lattice for different q < 0.5 
agree very well with this theory by reproducing the log- 
arithmic growth law (fH*)]) with its prefactor 2k [see fig. 
[5Jb)]. This is even true for q approaching 0.5 when the 
domain interface in fig. IDJd) becomes broader and, strictly 
speaking, the bounding wall of our model does not exhibit 
a hard-core repulsion for the random walker. Just as in 
the q = p case, simulations of an extending lattice are 
independent of a and 7. Note that the analogy with a 
biased random walker bounded by a wall also gives the 
growth laws N(t) oc t for q > p and N(t) cx \ft for q = p. 
We checked this by Monte-Carlo simulations. 

Summarizing our results, we found three distinct 
regimes for the bias q in the extending TASEP with qual- 
itatively different growth functions. The experimental ob- 
servation of logarithmic growth for flagellar filaments [5] 
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Fig. 7: Lattice growth speed from simulations for small neg- 
ative bias (q = 0.4999). The inset shows flagellar elongation 
data of three samples of Salmonella and fits, taken from [3]. 



mentioned in the introduction is in qualitative agreement 
with our model when the drift is negative (q < 0.5). In 
order to compare our result with the experiment in more 
detail, we carried out simulations with very small negative 
drift. By setting q = 0.4999, it was possible to achieve 
a lattice length of the order of N = 20000 in reason- 
able simulation time while still reaching the logarithmic 
growth regime. If each growth event mimics one flagcllin 
molecule, N = 20000 corresponds to a filament of about 
10^m. For our estimate, we used the length 55A for a 
flagcllin molecule and the fact that 11 flagcllin molecules 
make up the circular cross section of the hollow flagellar 
filament. The speed of elongation of flagellar filaments in 
the experiment was found to decrease exponentially with 
increasing length for filament lengths in the range from 
4 to 14/xm (see inset of fig. [7]) [3]- This is in agreement 
with our simulations as illustrated in fig. [71 However, 
it remains to establish an explanation why the flagellin 
in the hollow filament should exhibit a negative drift to- 
wards the cell body. It has been suggested that proton 
motive force (PMF) is biasing the Brownian motion of 
the flagellin translocation process [3S] . The details of this 
mechanism are, however, very unclear at this stage and 
should be subject to further research. 
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